Multimédia¶

Trabalho 1 - Compressão de imagem com o JPEG¶

FCTUC - LEI - 2021/2022¶

  • ##### João Silva - 2019217672
  • ##### Joel Oliveira - 2019227468
  • ##### Tomás Mendes - 2019232272
In [58]:
import numpy as np
from matplotlib import image
import matplotlib.pyplot as plt
import matplotlib.colors as clr
from scipy.fftpack import dct
from scipy.fftpack import idct
import pandas as pd
import PIL as p
from matplotlib.colors import Normalize
%matplotlib inline
In [59]:
"""
Returns an array of type:
 - [x, y] for grayscale images
 - [x, y, [R, G, B]] for RGB images
 - [x, y, [R, G, B, A]] for RGBA images
"""
def read_image(filename):
    img = image.imread(filename)
    return img

3.2 função para criar colormap entre duas cores¶

In [60]:
"""
Funcao para criar um colormap;
Parametros:
            1.cmin-> tuplo com valores (r,g,b) minimos
            2.cmax-> tuplo com valores (r,g,b) maximos
Devolve o objeto correspondente ao colormap
"""
def create_colormap(cmin : tuple[float], cmax: tuple[float]):
    return clr.LinearSegmentedColormap.from_list('', [cmin, cmax], N=256)

3.3 Função para visualizar a imagem com um certo colormap¶

In [61]:
"""
Função para mostrar uma imagem. Aceita um colormap definido pelo utilizador ou os do matplotlib
"""
    
def show_image(img, colormap = None):    
    plt.figure(figsize=(8,8))
    
    # Imagens com apenas uma coponenete: R, G, B ou Grayscale
    if len(img.shape) == 2:
        plt.imshow(img, cmap = colormap)
    else:
        if colormap != None:
            new_img = img[:, :, 0]
            plt.imshow(new_img, cmap = colormap)
        else:
            plt.imshow(img)

    plt.axis('off')
    plt.show()

Semana 1¶

Na diretória imagens estão presentes as imagens jpeg com baixa, média e alta qualidade, de acordo com o nome das pastas

  • Análise da compressão das várias imagens em jpeg
  • Converter a imagem de RGB para YCbCr
  • Padding da imagen para ficar um múltiplo de 16x16

3.1 Leitura da imagem¶

In [62]:
#img = read_image('imagens/peppers.bmp')
img = read_image('imagens/barn_mountains.bmp')
#img = read_image('imagens/logo.bmp')

show_image(img)

Tamanho dos ficheiros .bmp¶

Barn Peppers Logo
size 356.5KB 589.9KB 421.6KB


Foi utilizado o Adobe Photoshop 2019 para comprimir as imagens.

Tamanho dos ficheiros após compressão para JPEG¶

Barn Peppers Logo
Low 43.4KB 35.2KB 21.9KB
Medium 51.5KB 41.3KB 23.1KB
High 67.5KB 57.7KB 27.3KB


Rácio de compressão¶

Barn Peppers Logo
Low 8.2:1 16.7:1 19.3:1
Medium 6.9:1 14.3:1 18.3:1
High 5.2:1 10.2:1 15.4:1


Diferenças subjetivas jpeg/bmp¶

Barn Peppers Logo
Low não muito evidente muito evidente evidente
Medium nada evidente não muito evidente evidente
High nada evidente nada evidente nada evidente

Comparação dos resultados¶

A compressão é maior em imagens cujos pixels adjacentes têm uma cor semelhante.
Isto sugere que o codec do Jpeg utiliza compressão baseada em adjacência.

Visualmente, notam-se maiores diferenças em locais onde há alterações bruscas na cor do pixel.

In [63]:
low = image.imread('imagens/Low/peppers.jpg')
low.shape
Out[63]:
(384, 512, 3)
In [64]:
medium = image.imread('imagens/Medium/peppers.jpg')
medium.shape
Out[64]:
(384, 512, 3)
In [65]:
# ColorMaps
cm_gray = clr.LinearSegmentedColormap.from_list('gray', [(0,0,0), (1, 1, 1)], N = 256)
cm_red = clr.LinearSegmentedColormap.from_list('red', [(0,0,0), (1, 0, 0)], N = 256)
cm_green = clr.LinearSegmentedColormap.from_list('green', [(0,0,0), (0, 1, 0)], N = 256)
cm_blue = clr.LinearSegmentedColormap.from_list('blue', [(0,0,0), (0, 0, 1)], N = 256)


# Adaptado de https://jakevdp.github.io/PythonDataScienceHandbook/04.07-customizing-colorbars.html
cmap = plt.cm.get_cmap(cm_gray)
gray = cmap(np.arange(cmap.N))

cmap = plt.cm.get_cmap(cm_red)
red = cmap(np.arange(cmap.N))

cmap = plt.cm.get_cmap(cm_green)
green = cmap(np.arange(cmap.N))

cmap = plt.cm.get_cmap(cm_blue)
blue = cmap(np.arange(cmap.N))

fig, ax = plt.subplots(4, figsize=(10, 5),subplot_kw=dict(xticks=[], yticks=[]))
ax[0].imshow([gray], extent=[0, 10, 0, 1])
ax[1].imshow([red], extent=[0, 10, 0, 1])
ax[2].imshow([green], extent=[0, 10, 0, 1])
ax[3].imshow([blue], extent=[0, 10, 0, 1])
Out[65]:
<matplotlib.image.AxesImage at 0x7f522e8b61f0>

3.4 Função para separar imagem nos seus componentes R, G, B¶

In [66]:
"""
Separar uma imagem RGB nos seus componentes
"""
def separate_components(img):
    r = img[:, :, 0]
    g = img[:, :, 1]
    b = img[:, :, 2]
    
    return r, g, b

3.5 Visualização da imagem nos seus componentes R, G, B¶

In [67]:
r, g, b = separate_components(img)
show_image(r, cm_red)
show_image(g, cm_green)
show_image(b, cm_blue)
In [68]:
"""
Juntar as coponentes R, G e B para formar uma imagem
"""
def join_components(r, g, b):
    return np.dstack((r, g, b))
In [69]:
img = join_components(r, g, b)
show_image(img)

4.1 Função para padding da imagem¶

In [70]:
"""
Recebe uma imagem e altera as suas dimensões (m,n) para (16*p, 16*q).
Isto é realizado através da cópia da ultima coluna/linha até atingir o valor multiplo de 16.
Conta o numero de linhas/colunas adicionadas, (x,y) : 0<=x,y<=15.
Devolve um tuplo com as dimensoes incrementadas (x,y) e a imagem com as novas dimensoes (16*p, 16*q)
"""
def padding(img : np.array):
    img = img.copy()
    shape = img.shape
    
    x,y = (16-img.shape[0]%16)%16, (16-img.shape[1]%16)%16
    h_padding = np.repeat(img[-1:,:,:], x, axis = 0)
    img = np.concatenate((img, h_padding), axis = 0)

    v_padding = np.repeat(img[:,-1:,:], y, axis = 1)
    img = np.concatenate((img, v_padding), axis = 1)
    return shape, img
In [71]:
"""
Recebe uma imagem e as dimensões originais dela.
Faz slice da imagem, elimando todos os elementos incrementados no padding;
Devolve a imagem original
"""
def unpad(img : np.array, shape):
    img = img.copy()
    x,y,z = shape
    img = img[:x, :y]
    return img
In [72]:
show_image(img)
shape, pad_img = padding(img)
show_image(pad_img)
unpad_img = unpad(img, shape)
show_image(unpad_img)

# Verificação que fica igual ao original
print(np.all(np.equal(img, unpad_img)))

print(f"Original shape: {img.shape}")
print(f"Padding shape: {pad_img.shape}")
print(f"After Padding and Unpadding shape: {unpad_img.shape}")

True
Original shape: (297, 400, 3)
Padding shape: (304, 400, 3)
After Padding and Unpadding shape: (297, 400, 3)

5.1 Função de conversão RGB para YCbCr¶

In [73]:
"""
Converte imagem no formato RGB para imagem no formato yCbCr;
"""
def rgb_to_ycbcr(img : np.array):
    img.copy().astype(float)
    
    y_cb_cr_mat = np.array([ [0.299    , 0.587    , 0.114    ]
                            ,[-0.168736, -0.331264, 0.5      ]
                            ,[0.5      , -0.418688, -0.081312] ], dtype=float)
    
    y  = y_cb_cr_mat[0,0] * img[:,:,0] + y_cb_cr_mat[0,1] * img[:,:,1] + y_cb_cr_mat[0,2]*img[:,:,2]
    cb = y_cb_cr_mat[1,0] * img[:,:,0] + y_cb_cr_mat[1,1] * img[:,:,1] + y_cb_cr_mat[1,2]*img[:,:,2] + 128
    cr = y_cb_cr_mat[2,0] * img[:,:,0] + y_cb_cr_mat[2,1] * img[:,:,1] + y_cb_cr_mat[2,2]*img[:,:,2] + 128
    
    y_cb_cr = np.dstack((y, cb, cr))
    return y_cb_cr
In [74]:
"""
Converte imagem no formato YCbCr para imagem no formato RGB
"""
def ycbcr_to_rgb(img : np.array):
    img = img.copy().astype(float)
    
    y_cb_cr_mat_inv = np.linalg.inv(
                                np.array([ [0.299    , 0.587    , 0.114    ]
                                        ,  [-0.168736, -0.331264, 0.5      ]
                                        ,  [0.5      , -0.418688, -0.081312] ], dtype=float)
                                    )
    y = img[:, :, 0]
    cb = img[:, :, 1] - 128
    cr = img[:, :, 2] - 128
    
    r = y + y_cb_cr_mat_inv[0,2]*cr
    g = y + y_cb_cr_mat_inv[1,1]*cb + y_cb_cr_mat_inv[1,2]*cr
    b = y + y_cb_cr_mat_inv[2,1]*cb
    
    rgb = np.dstack((r,g,b))
    rgb = np.round(rgb)
    rgb[rgb > 255] = 255
    rgb[rgb < 0] = 0
    
    return rgb.astype(np.uint8)
In [75]:
ycbcr = rgb_to_ycbcr(pad_img)
rgb = ycbcr_to_rgb(ycbcr)

# Verificação que fica igual ao original
#print(np.all(np.equal(img, rgb)))

show_image(img)
show_image(rgb)
In [76]:
# Método do colormap para representar os canais Cb e Cr:
# https://stackoverflow.com/questions/28638848/displaying-y-cb-and-cr-components-in-matlab

cm_cb = clr.LinearSegmentedColormap.from_list('cb', [(0.5, 0.5, 0), (0.5, 0.5, 1)], N = 256)
cm_cr = clr.LinearSegmentedColormap.from_list('cr', [(0, 0.5, 0.5), (1, 0.5, 0.5)], N = 256)

cmap = plt.cm.get_cmap(cm_cb)
cm_cb_rep = cmap(np.arange(cmap.N))

cmap = plt.cm.get_cmap(cm_cr)
cm_cr_rep = cmap(np.arange(cmap.N))

fig, ax = plt.subplots(2, figsize=(10, 2),subplot_kw=dict(xticks=[], yticks=[]))
ax[0].imshow([cm_cb_rep], extent=[0, 10, 0, 1])
ax[1].imshow([cm_cr_rep], extent=[0, 10, 0, 1])
Out[76]:
<matplotlib.image.AxesImage at 0x7f522e9e82b0>
In [77]:
y, cb, cr = separate_components(ycbcr)

show_image(y, cm_gray)
show_image(cb, cm_cb)
show_image(cr, cm_cr)
In [78]:
f, ax = plt.subplots(nrows=3,ncols=2,figsize=(12,12))

ax[0,0].imshow(y, cm_gray)
ax[0,1].imshow(r, cm_red)

ax[1,0].imshow(cb, cm_cb)
ax[1,1].imshow(g, cm_green)

ax[2,0].imshow(cr, cm_cr)
ax[2,1].imshow(b, cm_blue)
for i in ax.flatten():
    i.set_axis_off()

Verifica-se que o canal Y da imagem tem muito mais definição do que o Cb e do que o Cr.

Em relação aos canais R,G e B, os canais R e G são os que apresentam maior semelhança visual com o canal Y. O canal G é o que tem mais peso no calculo do canal Y, seguido do R e finalmente do B.

Isto deve-se ao facto dos niveis de luminância ser maior precisamente nos canais G, seguido do R, seguido do B.
Com isto, concentra-se a maior parte da informação da luminância num dos canais, e informação sobre a crominância, que é menos importante para o olho humano nos outros (Cb e Cr).

Isto permite que haja maior compressão nos canais Cb e Cr, pois esta informação afeta menos a perceção do olho humano.

Semana 2¶

Aplicação dos primeiros passos destrutivos do CODEC jpeg

  • Downsample com os rácios 4:2:2 e 4:2:0
  • DCT aplicada á imagem toda

6.1 Downsample¶

In [79]:
"""
(Y:Cb:Cr)
(4, 2, 0) -> São removidos pixeis nas linhas e colunas
(4, 2, 2) -> São removidos pixeis apenas nas colunas

"""
def sub_sample(y, cb, cr, downsample_ratio):
    y = y.copy()
    cb = cb.copy().astype(float)
    cr = cr.copy().astype(float)
    if downsample_ratio == (4,2,0):
        ratio = round(downsample_ratio[0]/downsample_ratio[1])
        cb = cb[::ratio,::ratio]
        cr = cr[::ratio,::ratio]
    else:
        v_ratio = round(downsample_ratio[0]/downsample_ratio[1])
        h_ratio = round(downsample_ratio[0]/downsample_ratio[2])
        
        cb = cb[:, ::v_ratio]
        cr = cr[:, ::h_ratio]
    return y,cb,cr
In [80]:
"""
funciona em -- (4,2,0) e (4,2,2) --
Parametros:
            1. y --------------> matriz com os valores da luminancia Y
            2. cb -------------> matriz downsampled de cb
            3. cr -------------> matriz downsample de cr
            4. upsample_ratio -> tuplo com 3 inteiros correspondentes ao ratio de downsample (Y:Cb:Cr)
"""
def const_up_sample(y, cb, cr, upsample_ratio):
    y = y.copy()
    cb = cb.copy()
    cr = cr.copy()
    
    cbArr = np.zeros(shape = y.shape, dtype=float)
    crArr = np.zeros(shape = y.shape, dtype=float)
    
    if upsample_ratio[-1] == 0:
        ratio = int(upsample_ratio[0]/upsample_ratio[1])
        
        cbArr[::ratio,::ratio] = cb
        cbArr[1::ratio,::ratio] = cb
        cbArr[::ratio,1::ratio] = cb
        cbArr[1::ratio,1::ratio] = cb
        
        crArr[::ratio,::ratio] = cr
        crArr[1::ratio,::ratio] = cr
        crArr[::ratio, 1::ratio] = cr
        crArr[1::ratio,1::ratio] = cr
        
    else:
        cb_ratio = int(upsample_ratio[0]/upsample_ratio[1])
        cr_ratio = int(upsample_ratio[0]/upsample_ratio[2])
        
        cbArr[:,::cb_ratio] = cb
        cbArr[:,1::cb_ratio] = cb
        
        crArr[:,::cr_ratio] = cr
        crArr[:,1::cr_ratio] = cr
        
    return y,cbArr,crArr
In [81]:
"""
funciona em -- (4,2,0) e (4,2,2) --
Parametros:
            1. y --------------> matriz com os valores da luminancia Y
            2. cb -------------> matriz downsampled de cb
            3. cr -------------> matriz downsample de cr
            4. upsample_ratio -> tuplo com 3 inteiros correspondentes ao ratio de downsample (Y:Cb:Cr)
"""
def linear_up_sample(y, cb, cr, upsample_ratio):
    y = y.copy()
    cb = cb.copy()
    cr = cr.copy()
    
    new_cb = np.zeros(shape = y.shape, dtype=float)
    new_cr = np.zeros(shape = y.shape, dtype=float)
    
    if upsample_ratio[2] == 0:
        ratio = int(upsample_ratio[0]/upsample_ratio[1])
        #--------------------Cb interpolation----
        cb_cols_mean = (cb[:,:-1] + np.roll(cb, -1, 1)[:,:-1])//2
        cb_cols_mean = np.concatenate((cb_cols_mean, cb[:,-1:]), axis = 1)

        new_cb[::2,::2] = cb
        new_cb[::2,1::2] = cb_cols_mean

        cb = new_cb[::2]
        
        cb_lines_mean = (cb[:-1] + np.roll(cb, -1, 0)[:-1])//2
        cb_lines_mean = np.concatenate((cb_lines_mean, cb[-1:]), axis = 0)

        new_cb[1::2,:] = cb_lines_mean
        
        #------------------------Cr interpolation----
        cr_cols_mean = (cr[:,:-1] + np.roll(cr, -1, 1)[:,:-1])//2
        cr_cols_mean = np.concatenate((cr_cols_mean, cr[:,-1:]), axis = 1)

        new_cr[::2,::2] = cr
        new_cr[::2,1::2] = cr_cols_mean

        cr = new_cr[::2]
        
        cr_lines_mean = (cr[:-1] + np.roll(cr, -1, 0)[:-1])//2
        cr_lines_mean = np.concatenate((cr_lines_mean, cr[-1:]), axis = 0)

        new_cr[1::2,:] = cr_lines_mean
        
    else:
        cb_ratio = int(upsample_ratio[0]/upsample_ratio[1])
        cr_ratio = int(upsample_ratio[0]/upsample_ratio[2])
        
        cb_mean = (cb[:,:-1] + np.roll(cb, -1, 1)[:,:-1])/2
        cb_mean = np.concatenate((cb_mean, cb[:,-1:]), axis = 1)
        new_cb[:,::cb_ratio] = cb
        new_cb[:,1::cb_ratio] = cb_mean
        
        cr_mean = (cr[:,:-1] + np.roll(cr, -1, 1)[:,:-1])/2
        cr_mean = np.concatenate((cr_mean, cr[:,-1:]), axis = 1)
        new_cr[:,::cr_ratio] = cr
        new_cr[:,1::cr_ratio] = cr_mean
    
    return y,new_cb, new_cr
In [82]:
ratio = (4, 2, 2)

y_d, cb_d, cr_d = sub_sample(y, cb, cr, ratio)

show_image(y_d, cm_gray)
show_image(cb_d, cm_cb)
show_image(cr_d, cm_cr)

print("Y shape: " + str(y_d.shape), 
      "Cb shape: " + str(cb_d.shape), 
      "Cr shape: " + str(cr_d.shape), 
      sep = '\n'
     )

Y shape: (304, 400)
Cb shape: (304, 200)
Cr shape: (304, 200)
In [83]:
ratio = (4, 2, 0)

y_d, cb_d, cr_d = sub_sample(y, cb, cr, ratio)

show_image(y_d, cm_gray)
show_image(cb_d, cm_cb)
show_image(cr_d, cm_cr)

print("Y shape: " + str(y_d.shape), 
      "Cb shape: " + str(cb_d.shape), 
      "Cr shape: " + str(cr_d.shape), 
      sep = '\n'
     )

Y shape: (304, 400)
Cb shape: (152, 200)
Cr shape: (152, 200)

6.3 Análise da aplicação do downsample¶

No downsampling 4:2:2, há uma redução do número de pixeis, nos canais Cb e Cr, de 2:1. No caso do 4:2:0, a taxa de compressão é maior, pois é de 4:1 nesses mesmos canais.

Esta etapa de downsampling é destrutiva, pelo que a contrapartida é que quanto mais alta for a compressão, mais alta é a taxa de destrutividade, pois quantos mais dados forem eliminados, menos fiel vai ser a reconstrução (tendo em conta os nosso preditores).

Assim, a escolha destas métricas deve ter em conta a qualidade que se pretende obter na reconstrução da imagem.

7.1.1 DCT transformation¶

In [84]:
def dct_(a):
    return dct(dct(a, norm="ortho").T, norm = "ortho").T

7.1.1 DCT reverse transformation¶

In [85]:
def idct_(a_dct):
    return idct(idct(a_dct, norm='ortho').T, norm='ortho').T
In [86]:
# DCT aplicada na imagem toda

y_dct = dct_(y_d)
cb_dct = dct_(cb_d)
cr_dct = dct_(cr_d)

show_image(np.log(np.abs(y_dct) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dct) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dct) + 0.0001), cm_gray)
In [87]:
# Inverso da DCT

y_idct = idct_(y_dct)
cb_idct = idct_(cb_dct)
cr_idct = idct_(cr_dct)

show_image(y_idct, cm_gray)
show_image(cb_idct, cm_gray)
show_image(cr_idct, cm_gray)
In [88]:
# Reconstrução com copia da coluna da esquerda

y, cb, cr = const_up_sample(y_idct, cb_idct, cr_idct, ratio)

show_image(y, cm_gray)
show_image(cb, cm_gray)
show_image(cr, cm_gray)
In [89]:
#Recontrução com interpolação linear (média) da coluna anterior e seguinte

y, cb, cr = linear_up_sample(y_idct, cb_idct, cr_idct, ratio)

show_image(y, cm_gray)
show_image(cb, cm_gray)
show_image(cr, cm_gray)

7.1.3 Potencial de Compressão¶

A DCT, tal como a DFT, permite representar um sinal por valores de frequência com coeficientes associados. No entanto, a DCT tem a vantagem de conseguir concentrar muito mais a energia.

No caso, a DCT faz uma transformação do domínio do espaço para o domínio da frequência. Nesta transformação, os coeficientes das baixas frequências ficam presentes no canto superior esquerdo e os das altas frequências, no canto inferior direito.

Normalmente, no amplo espaço de uma imagem há vários tipos de transições, abruscas e suaves. Como é visivel em todas as imagens não há uma concentração de energia notória na DCT das imagens. Aliás, a matriz resultante demonstra-se semelhante ao ruído branco, que tem uma baixa taxa de compressão.

Semana 3¶

Aplicação da DCT em blocos 8x8 ou 64x64

  • DCT em blocos e inverso

7.2.1 Dct in blocks¶

In [90]:
def dct_n_blocks(y, cb, cr, step):
    dct_y = np.zeros(y.shape)
    dct_cb = np.zeros(cb.shape)
    dct_cr = np.zeros(cr.shape)
    
    for i in range(0, y.shape[0], step):
        for j in range(0, y.shape[1], step):
            dct_y[i:i+step, j:j+step] = dct_(y[i:i+step, j:j+step])
    
    for i in range(0, cb.shape[0], step):
        for j in range(0, cb.shape[1], step):
            dct_cb[i:i+step, j:j+step] = dct_(cb[i:i+step, j:j+step])
            dct_cr[i:i+step, j:j+step] = dct_(cr[i:i+step, j:j+step])
            
    return dct_y, dct_cb, dct_cr

7.2.1 Reverse DCT in blocks¶

In [91]:
def idct_n_blocks(y, cb, cr, step):
    idct_y = np.zeros(y.shape)
    idct_cb = np.zeros(cb.shape)
    idct_cr = np.zeros(cr.shape)
    
    for i in range(0, y.shape[0], step):
        for j in range(0, y.shape[1], step):
            idct_y[ i : i+step, j:j+step] = idct_(y[i:i+step, j:j+step])
    
    for i in range(0, cb.shape[0], step):
        for j in range(0, cb.shape[1], step):
            idct_cb[ i : i+step, j : j+step] = idct_(cb[ i : i+step, j : j+step])
            idct_cr[ i : i+step, j : j+step] = idct_(cr[ i : i+step, j : j+step])
            
    return idct_y, idct_cb, idct_cr
In [92]:
# DCT em blocos

n = 64
y_dctn, cb_dctn, cr_dctn = dct_n_blocks(y, cb, cr, n)

show_image(np.log(np.abs(y_dctn) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn) + 0.0001), cm_gray)
In [93]:
# Inverter DCT em blocos

y_d, cb_d, cr_d = idct_n_blocks(y_dctn, cb_dctn, cr_dctn, n)

show_image(y_d, cm_gray)
show_image(cb_d, cm_gray)
show_image(cr_d, cm_gray)
In [94]:
# DCT em blocos

n = 8
y_dctn, cb_dctn, cr_dctn = dct_n_blocks(y, cb, cr, n)

show_image(np.log(np.abs(y_dctn) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn) + 0.0001), cm_gray)
In [95]:
# Inverter DCT em blocos

y_d, cb_d, cr_d = idct_n_blocks(y_dctn, cb_dctn, cr_dctn, n)

show_image(y_d, cm_gray)
show_image(cb_d, cm_gray)
show_image(cr_d, cm_gray)

7.2.3/7.3.2 DCT aplicada a blocos 8x8 e 64x64¶

A aplicação da DCT por bloco parte do princípio que num pequeno espaço da imagem, a variação vai ser muito baixa.

Ao contrário do que aconteceu no cenário anterior, em que foi aplicada a DCT em toda a imagem, no caso da aplicação da DCT em blocos 8x8 é clara a concentração da energia em cada bloco, o que aumenta o potencial de compressão.
No caso da transformada aplicada a blocos 64x64, também é visível a concentração da energia, mas não tanto como no caso 8x8. É visível, em cada bloco, algo semelhante a ruído branco.

Semana 4¶

Quantização da imagem e compressão dos coponentes DC da DCT com delta encoding

  • Quantização dos canais Y, Cb e Cr com as matrizes adequadas
  • delta encoding para as componentes DC da DCT.
    • A compressão dos componentes AC não faz parte deste trabalho

8 Quantização¶

In [96]:
def get_Qs():
    Q_y = np.array(
        [
            [16, 11, 10, 16, 24, 40, 51, 61],
            [12, 12, 14, 19, 26, 58, 60, 55],
            [14, 13, 16, 24, 40, 57, 69, 56],
            [14, 17, 22, 29, 51, 87, 80, 62],
            [18, 22, 37, 56, 68, 109, 103, 77],
            [24, 35, 55, 64, 81, 104, 113, 92],
            [49, 64, 78, 87, 103, 121, 120, 101],
            [72, 92, 95, 98, 112, 100, 103, 99],
        ]
    )
    
    Q_CbCr = np.array(
        [
            [17, 18, 24, 47, 99, 99, 99, 99],
            [18, 21, 26, 66, 99, 99, 99, 99],
            [24, 26, 56, 99, 99, 99, 99, 99],
            [47, 66, 99, 99, 99, 99, 99, 99],
            [99, 99, 99, 99, 99, 99, 99, 99],
            [99, 99, 99, 99, 99, 99, 99, 99],
            [99, 99, 99, 99, 99, 99, 99, 99],
            [99, 99, 99, 99, 99, 99, 99, 99],
        ]
    )
    return Q_y, Q_CbCr
In [97]:
def quantitization(y, cb, cr, fator):
    q_y, q_cbcr = get_Qs()
    y_q = np.zeros(y.shape)
    cb_q = np.zeros(cb.shape)
    cr_q = np.zeros(cr.shape)
    
    if fator>=50:
        s = (100-fator)/50
    else:
        s = 50/fator
    
    if s!=0:
        q_y = np.round(s*q_y)
        q_cbcr = np.round(s*q_cbcr)
    else:
        q_y[::,::] = 1
        q_cbcr[::,::] = 1
    
    q_y[q_y>255] = 255    
    q_cbcr[q_cbcr>255] = 255
    
    for r in range(0,y.shape[0],8):
        for c in range(0,y.shape[1],8):
            y_q[r:r+8, c:c+8] = y[ r : r+8, c : c+8] / q_y
            
    for r in range(0, cb.shape[0], 8):
        for c in range(0,cb.shape[1], 8):
            cb_q[ r : r+8, c : c+8] = cb[ r : r+8, c : c+8] / q_cbcr
            cr_q[ r : r+8, c : c+8] = cr[ r : r+8, c : c+8] / q_cbcr
    
    y_q = np.round(y_q)
    cb_q = np.round(cb_q)
    cr_q = np.round(cr_q)
    return y_q, cb_q, cr_q
In [98]:
def quantitization_inverse(y, cb, cr, fator):
    q_y, q_cbcr = get_Qs()
    y_q = np.zeros(y.shape, dtype=float)
    cb_q = np.zeros(cb.shape, dtype= float)
    cr_q = np.zeros(cr.shape, dtype=float)
    
    if fator>=50:
        s = (100-fator)/50
    else:
        s = 50/fator
    
    if s!=0:
        q_y = np.round(s*q_y)
        q_cbcr = np.round(s*q_cbcr)
    else:
        q_y[::,::] = 1
        q_cbcr[::,::] = 1
    
    q_y[ q_y > 255 ] = 255    
    q_cbcr[ q_cbcr > 255 ] = 255
    
    for r in range(0, y.shape[0], 8):
        for c in range(0, y.shape[1], 8):
            y_q[ r : r+8, c : c+8] = y[ r : r+8, c : c+8]*q_y
            
    for r in range(0, cb.shape[0], 8):
        for c in range(0, cb.shape[1], 8):
            cb_q[ r : r+8, c : c+8] = cb[ r : r+8, c : c+8] * q_cbcr
            cr_q[ r : r+8, c : c+8] = cr[ r : r+8, c : c+8] * q_cbcr
    
    y_q = np.round(y_q)
    cb_q = np.round(cb_q)
    cr_q = np.round(cr_q)
    return y_q, cb_q, cr_q
In [99]:
# Fator de qualidade para a quantização da imagem

fator = 100
y_dctn_q, cb_dctn_q, cr_dctn_q = quantitization(y_dctn, cb_dctn, cr_dctn, fator)
y_dctn_qi, cb_dctn_qi, cr_dctn_qi = quantitization_inverse(y_dctn_q, cb_dctn_q, cr_dctn_q, fator)

show_image(np.log(np.abs(y_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn_q) + 0.0001), cm_gray)
In [100]:
# Fator de qualidade para a quantização da imagem

fator = 75
y_dctn_q, cb_dctn_q, cr_dctn_q = quantitization(y_dctn, cb_dctn, cr_dctn, fator)
y_dctn_qi, cb_dctn_qi, cr_dctn_qi = quantitization_inverse(y_dctn_q, cb_dctn_q, cr_dctn_q, fator)

show_image(np.log(np.abs(y_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn_q) + 0.0001), cm_gray)
In [101]:
# Fator de qualidade para a quantização da imagem

fator = 50
y_dctn_q, cb_dctn_q, cr_dctn_q = quantitization(y_dctn, cb_dctn, cr_dctn, fator)
y_dctn_qi, cb_dctn_qi, cr_dctn_qi = quantitization_inverse(y_dctn_q, cb_dctn_q, cr_dctn_q, fator)

show_image(np.log(np.abs(y_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn_q) + 0.0001), cm_gray)
In [102]:
# Fator de qualidade para a quantização da imagem

fator = 25
y_dctn_q, cb_dctn_q, cr_dctn_q = quantitization(y_dctn, cb_dctn, cr_dctn, fator)
y_dctn_qi, cb_dctn_qi, cr_dctn_qi = quantitization_inverse(y_dctn_q, cb_dctn_q, cr_dctn_q, fator)

show_image(np.log(np.abs(y_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn_q) + 0.0001), cm_gray)
In [103]:
# Fator de qualidade para a quantização da imagem

fator = 10
y_dctn_q, cb_dctn_q, cr_dctn_q = quantitization(y_dctn, cb_dctn, cr_dctn, fator)
y_dctn_qi, cb_dctn_qi, cr_dctn_qi = quantitization_inverse(y_dctn_q, cb_dctn_q, cr_dctn_q, fator)

show_image(np.log(np.abs(y_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn_q) + 0.0001), cm_gray)

8.3 Potencial de compressão em função do fator de qualidade¶

É possível verificar que à medida que o fator de qualidade diminuiu, também diminuiu o número de frequências que reprentam cada bloco. Em casos em que o fator é reduzido há um grande número de blocos que aparentam ser representados apenas pelo seu componente DC.

Com a imagem representada por um menor número de diferentes valores, o potencial de compressão aumenta, pois a distribuição dos valores diminuiu, tal como a entropia. Assim, quanto menor for o fator de qualidade, maior é o potencial de compressão.

8.4 Comparação com resultados da alínea 7¶

Na alínea 5 foi referido que o canal Y é calculado de forma a ter a maior parte da informação de luminância da imagem e que os canais Cb e Cr são calculados de modo a conterem mais informação de crominância.
Expôs-se também que a crominância tem menos relevância na perceção da imagem do olho humano, e, portanto, os canais Cb e Cr teriam mais liberdade para a perda de informação.

É possível notar-se em todos os exemplos que o canal Y mantem sempre mais informação que os canais Cb e Cr. Isto deve-se exatamente à tentativa de minimizar a perda de informação de luminância na compressão destrutiva.

9.2 DPCM¶

In [104]:
# Delta encoding aplicado aos coeficientes DC da DCT em blocos

def diff(bloco):
    x, y = bloco.shape
    bloco = bloco.flatten().astype(float)
    bloco[1:] -= bloco[:-1]
    return bloco.reshape(x,y)

def DPCM(y,cb,cr):

    yBlocks = y.copy()
    cbBlocks = cb.copy()
    crBlocks = cr.copy()
    
    yBlocks[::8,::8] = diff(y[::8,::8])
    cbBlocks[::8,::8] = diff(cb[::8,::8])
    crBlocks[::8,::8] = diff(cr[::8,::8])
        
    return yBlocks, cbBlocks, crBlocks
In [105]:
def diff_reverse(bloco):
    x, y = bloco.shape
    bloco = bloco.flatten().astype(float)
    for i in range(1,bloco.shape[0]):
        bloco[i] = bloco[i] + bloco[i-1]
    return bloco.reshape(x,y)

def DPCM_reverse(yBlocks, cbBlocks, crBlocks):
    y = yBlocks.copy()
    cr = crBlocks.copy()
    cb = cbBlocks.copy()
    
    y[::8,::8] = diff_reverse(yBlocks[::8,::8])
    cr[::8,::8] = diff_reverse(crBlocks[::8,::8])
    cb[::8,::8] = diff_reverse(cbBlocks[::8,::8])
    
    return y, cb, cr
In [106]:
d_y, d_cb, d_cr = DPCM(y_dctn_q, cb_dctn_q, cr_dctn_q)
id_y, id_cb, id_cr = DPCM_reverse(d_y, d_cb, d_cr)

show_image(np.log(np.abs(d_y) + 0.0001), cm_gray)
show_image(np.log(np.abs(d_cb) + 0.0001), cm_gray)
show_image(np.log(np.abs(d_cr) + 0.0001), cm_gray)

Análise sobre a aplicação da codificação DPCM¶

Após a aplicação da codificação DPCM, o número de blocos contíguos totalmente nulos aumentou bastante. Isto potencializa a compressão de algorítmos como o RLE.

De notar que esta compressão não é destrutiva, pelo que este passo ajuda a compressão sem perder qualidade na reconstrução da imagem, através de um método simples de aplicar.

Semana 5¶

Aplicação de todas as funções criadas antriormente com os fatores de qualidade 10, 25, 50, 75, 100

  • Cálculo das métricas de distorção (MSE, RMSE, SNR, PSNR) para todas as imagens e fatores de qualidade
In [107]:
def MSE(orig, compressed):
    return np.sum(np.square(orig-compressed))/np.prod(orig.shape) * 3
In [108]:
def RMSE(orig, compressed):
    return np.sqrt(MSE(orig,compressed))
In [109]:
def SNR(orig, compressed):
    P = np.sum(np.square(orig))/np.prod(orig.shape)
    return 10*np.log10(P/MSE(orig, compressed))
In [110]:
def PSNR(orig, compressed):
    return 10*np.log10(np.square(orig).max()/MSE(orig,compressed))

Encoder and decoder functions¶

In [111]:
def encoder(filename, ratio, n, factor):
    img = read_image(filename)
    
    shape, pad_img = padding(img)
    
    ycbcr = rgb_to_ycbcr(pad_img)
    
    y, cb, cr = separate_components(ycbcr)
    
    y_d, cb_d, cr_d = sub_sample(y, cb, cr, ratio)
    
    y_dct, cb_dct, cr_dct = dct_n_blocks(y_d, cb_d, cr_d, n)
    
    y_dct_q, cb_dct_q, cr_dct_q = quantitization(y_dct, cb_dct, cr_dct, factor)
    
    delta_y, delta_cb, delta_cr = DPCM(y_dct_q, cb_dct_q, cr_dct_q)
    
    return delta_y, delta_cb, delta_cr, shape
In [112]:
def decoder(delta_y, delta_cb, delta_cr, ratio, shape, n, factor):
    y_dct_q, cb_dct_q, cr_dct_q = DPCM_reverse(delta_y, delta_cb, delta_cr)
    
    y_dct, cb_dct, cr_dct = quantitization_inverse(y_dct_q, cb_dct_q, cr_dct_q, factor)
    
    y_idct, cb_idct, cr_idct = idct_n_blocks(y_dct, cb_dct, cr_dct, n) # DCT em blocos n x n
    
    #y, cb, cr = const_up_sample(y_idct, cb_idct, cr_idct, ratio) # Upsample: cópia da coluna á esquerda
    y, cb, cr = linear_up_sample(y_idct, cb_idct, cr_idct, ratio) # Upsample: média entre o anterior e o seguinte
    
    ycbcr = join_components(y, cb, cr)
    
    img_rgb = ycbcr_to_rgb(ycbcr)
    
    img = unpad(img_rgb, shape)
    
    return img
In [113]:
%matplotlib inline

f = "imagens/barn_mountains.bmp"

ratio = (4,2,2)
n = 8
fator = 75
x,y,z,shape = encoder(f,ratio,n,fator)
img = decoder(x,y,z,ratio,shape,n,fator)
show_image(img)

print("MSE:", MSE(read_image(f).astype(float), img.astype(float)),
     "RMSE:", RMSE(read_image(f).astype(float), img.astype(float)),
     "SNR:", SNR(read_image(f).astype(float), img.astype(float)),
     "PSNR:",PSNR(read_image(f).astype(float), img.astype(float)),
     sep = "\n")

MSE:
152.02205387205387
RMSE:
12.329722376114308
SNR:
20.595997303892048
PSNR:
26.311737651588494
In [114]:
imagens = ['imagens/barn_mountains.bmp', 'imagens/peppers.bmp', 'imagens/logo.bmp']

ratio = (4,2,2)
n = 8
fatores = [10, 25, 50, 75, 100]

for file in imagens:
    print(f"{file.split('/')[1]}")
    f, ax = plt.subplots(nrows = 1, ncols = 5, figsize = (20, 5))
    f.suptitle(f"{file.split('/')[1]}", fontsize = 16)
    
    df = pd.DataFrame(index = fatores, columns = ["MSE", "RMSE", "SNR", "PSNR"])
    
    for index, fator in enumerate(fatores):
        y, cb, cr, shape = encoder(file, ratio, n, fator)
        i_jpeg = decoder(y, cb, cr, ratio, shape, n, fator)

        ycbcr = rgb_to_ycbcr(i_jpeg)
        y, cb, cr = separate_components(ycbcr)

        img = read_image(file)
        ycbcr = rgb_to_ycbcr(img)
        y_orig, cb_orig, cr_orig = separate_components(ycbcr)

        y_diff = np.abs(y_orig - y)
        ax[index].title.set_text(f"{file.split('/')[1]}: {fator}")
        y_diff[-1,-1] = 255
        y_diff[-1, -2] = 0 
        ax[index].imshow(y_diff, cm_gray)
        
        df["MSE"][fator] = MSE(img.astype(float), i_jpeg.astype(float))
        df["RMSE"][fator] =  RMSE(img.astype(float), i_jpeg.astype(float))
        df["SNR"][fator] = SNR(img.astype(float), i_jpeg.astype(float))
        df["PSNR"][fator] = PSNR(img.astype(float), i_jpeg.astype(float))
    
    print(df)
    for i in ax:
        i.set_axis_off()
    
    

barn_mountains.bmp
            MSE       RMSE        SNR       PSNR
10    707.49287  26.598738  13.917843  19.633583
25   398.153131  19.953775  16.414562  22.130302
50   261.369487  16.166926  18.242514  23.958255
75   152.022054  12.329722  20.595997  26.311738
100   11.170572   3.342241  31.934309  37.650049
peppers.bmp
            MSE       RMSE        SNR       PSNR
10   283.115982  16.826051  15.635096   23.61116
25   127.319875  11.283611  19.105778  27.081842
50    77.324849   8.793455  21.271549  29.247613
75     50.02124    7.07257  23.163195  31.139259
100    6.346832   2.519292  32.129169  40.105233
logo.bmp
            MSE       RMSE        SNR       PSNR
10   157.304726  12.542118  24.593488  26.163386
25    65.218448   8.075794  28.417201  29.987099
50    43.806448   6.618644  30.145526  31.715423
75    23.667495   4.864925  32.819383  34.389281
100    6.184128   2.486791  38.648121  40.218019

10.2 Métricas de distorção¶

Tanto o ruído quanto a diferença entre os pixeis da imagem original e da imagem reconstruida serão tanto maiores quanto menor for fator de qualidade. Desta forma, pode-se concluir que o SNR e o PSNR terão valores superiores em imagens cujo fator de qualidade seja alto e o MSE e RMSE terão valores tanto maiores quanto menor for o fator de qualidade.

In [ ]: